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^ Abstract 

o 

^ We provide a general method to compute a Taylor expansion in time of implied volatility for 

• stochastic volatility models, using a heat kernel expansion. Beyond the order implied volatility 

which is already known, wc compute the first order correction exactly at all strikes from the 

^ scalar coefficient of the heat kernel expansion. Furthermore, the first correction in the heat kernel 

expansion gives the second order correction for implied volatility, which we also give exactly at all 
strikes. As an application, we compute this asymptotic expansion at order 2 for the SABR model. 
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1 Introduction 

The most known model for pricing derivatives is the Black-Scholes-Merton model, where the under- 
lying is supposed to follow a geometric Brownian motion. Popular extensions include local volatility 
models and stochastic volatility models. For example the SABR model |HKLW02] combines the local 
volatility of the CEV model |Cox75] and a lognormal volatility process. Closed formulas for European 
options can be obtained for a few models; it is the case of the CEV model or for a stochastic volatility 
example the Heston model |Hes93j . These are however special cases and there are generally no closed 
form formulas. Finite difference methods or Monte-Carlo simulations can be used to price derivatives. 
Approximations has also been computed to achieve faster pricing, especially for calibration processes. 

For short maturities, Hagan, Kumar, Lesniewski and Woodward provide an approximation for 
the implied volatility of the SABR model they introduce [HKLW02, . Berestycki, Busca and Flo- 
rent |BBF02| IBBF04| and Henry-Labordere |HL05] give general methods to compute short maturity 
asymptotics of stochastic volatility models. These expansions give the implied volatility at first order 
in maturity. In addition some quantities are approximated by their value at the money, which can 
produce errors in the wings of the distributions. 

In this paper, we leverage on the heat kernel methods introduced for the study of stochastic 
processes by Varadhan |Var67al IVar67bj and used for the SABR model in |HLW05j and |HL05j . Using 
the heat kernel expansion of DeWitt |DeW65j . we provide a method to compute exactly a Taylor 
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expansion of the implied volatility at all strikes. The stochastic volatility diffusion is formulated 
as a diffusion on a Riemannian manifold. The geodesic distance gives the implied volatility at null 
maturity. The multiplicative factor of the heat kernel expansion provides the first order (in time) 
correction to implied volatility. The first corrective term of the heat kerned is translated into the 
second order correction to implied volatility and similarly for higher order corrections. We perform 
a detailed computation up to order 2 of the Taylor expansion in time of implied volatility, without 
other approximations. 

More generally, our method can be used to approximate a stochastic volatility model by an other 
model for which a closed form solution exists, with an implied parameter computed as a Taylor 
expansion. 

As an application, we compute the asymptotic SABR volatility at order 2 and compare it to finite 
difference method results and to the original SABR expansion. 

Our results can be useful for pricing short maturities options or even long maturities options with 
very low volatility of volatility. When the approximation is not valid, a numerical method such as 
a finite difference method (FDM) has to be used. When it is valid, it gives much faster results. At 
very short maturities, the prices are even more precise. Calibration at short maturities appears to be 
more stable using this approximation. 

In section [2] we recast the financial model in physical and geometric terms and fix our conventions. 
In section [3] we use a heat kernel expansion to compute a short maturity expansion of Black or more 
generally CEV implied volatility. Finally in section |4] we apply the method to the SABR model and 
compare the results to FDM and to the original formula. 

2 Diffusion equation in covariant form 

A stochastic volatility model for some asset with pure diffusion (no jumps) is described by two risk- 
neutral processes: the asset price S and a variable V which describes the stochastic part of volatility. 
In Heston model V would be the variance whereas in the SABR model it is a factor of volatility. The 
diffusion is given by the stochastic differential equations 



where dWi and dW2 are two standard Brownian processes with correlation p. The dependence of 
parameters in variables S and V we have written is the more common, it may be more general with 
all parameters depending on both variables. 

Stochastic volatility models can be seen as diffusion on a Riemann surface. More precisely, prices 
of securities are sections of a line bundle over this Riemann surface which are solutions of a diffusion 
(or heat) equation. 

A introduction to this subject and its applications to finance can be found in j HL08] . We present 
here the formalism and define all quantities we use in order to set our conventions. 

2.1 Diffusion equation 

Let us consider a general model with n state variables X^(t) (which will be the spot and the volatility) 
which follow a pure diffusion process, without jumps. For simplicity we consider a European payoff 
of some maturity T. The price F{X{t),t) of such a payoff is the solution of a diffusion equation 



dS = fis{S)dt + as{S,V)dWi 
dV = fiv{V)dt + av{V)dW2 





(2) 
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where S*-' is the covariance matrix, /ij the drifts and r the numeraire rate. Ah coefficients can depend 
on state variables X^(t) and time t. Unless explicitly staten, we take Einstein sum convention: 
repeated indices are summed. The price of the European option is given by the solution of this 
equation with terminal boundary condition at maturity T given by the payoff. 

The covariance matrix S*-' can be seen mathematically as the inverse g^^ = S*-' of a metric gij on 
the space of variables. The diffusion equation describes the diffusion over a Riemannian manifold: 
the state of variables endowed with the metric gij = 

Examples 

1. The Black-Scholes equation in the monetary account numeraire with volatility a, risk free rate 
r and dividend yield q reads 



-dtF = (r - q)SdsF + ^-a^S'^dsdsF - rF 



r = r. 



2 

This is equation Q with = [r - q)S, T.^^ = a^S"^ and 
2. For the stochastic volatility model described by equations ([T|, /x* is a two-dimensional vector 

fJ-s 



The covariance matrix S*-' is 



pa say (Jv"^ 



In what follows we restrict ourselves to the case of time- homogeneous models: there is no explicit 
time-dependence in parameters. The generalization to time-dependent cases is not difficult. 

2.2 Gauge structure 

There are several gauge transformations which are natural for such systems: 

1. Change of numeraire: 

^(^''^ ^ Hx:r) 

where ^{X, t) is the price of a security which is always nonzero. Mathematically it is a real 
function which is positive everywhere, that we denote thus by ^{X,t) = e~^^'^'^\ 

2. Change of variables 

X — y X'{X) . 

The natural way to handle a system with gauge freedom is to introduce covariant derivatives. The 
coordinate freedom is handled through the Levi-Civita connection which acts respectively on scalars, 
vectors and 1-forms as 

Dif = dif 
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where F^^ are the Christoffel symbols. The action on tensors with more indices is obtained by acting 
on ah indices with the Christoffel symbols. Christoffel symbols can be computed from the metric as 

^ij = {di9ij + djgu - diQij) . 
A fundamental property of the Levi-Civita connection is the covariance of the metric: 

DiQjk = . 

The metric is used to transforms vectors into 1-forms and conversely, i.e. lowering or raising indices: 

A' = g^'A, 
At = ffijA'^ . 

The numeraire gauge freedom is handled through a M- valued connection!^ with spatial and time 
components given by a 1-form Ai and a scalaij^Q: 

ViF = {Di-Ai)F 
VtF = {dt-Q)F. 



Under the change of numeraire 



these operators are covariant, 



ViF — > e'^(^'*)ViF 
VtF — > e*(^'*)VjF , 



provided that Ai and Q are shifted as 



Ai — > Ai - di(t) 
Q Q-dt(P. 



Using these connections, the diffusion equation ([2| can be rewritten as 



-VtF = ^V'V^F . (3) 



Identifying terms between equations ([2]) and ([S]), the M connection must be 



A^ = <7^,(-^ri/'-^^) (4) 
Q = -^gi^(^diAj-AiAj-T'ljAk)+r. (5) 



In addition with 

gij = j^ij 

this translates the set of financial parameters into mathematical quantities. 



^ This connection is similar to tlie connection which described the electromagnetic potential, except that the fibre 
of the gauge bundle is R instead of (7(1). This causes a difference of a factor i in equations. 

^ There is a breaking of symmetry between time and spatial directions. The diffusion equation can be seen as a 
non-relativistic limit of a pure wave equation. 
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2.3 Kolmogorov forward equation 

The Kolmogorov backward equation ([S]) leads to a dual Kolmogorov forward equation. 

We suppose that all prices are expressed with respect to a numeraire which is a traded asset that 
does not pay any coupon or dividend. The price of the numeraire security itself is identically 1; this 
reads mathematically 

-Vti = ^v^va 

If p{X, t) is a risk-neutral probability density to get in state X at time t starting from state Xq 
at time 0, then the price of a European payoff of maturity T >t can be written as 

F(Xo,0) = j dXp{X,t)F{X,t) 

As t does not appear on the left-hand side, the derivative of the integral with respect to t must vanish. 



J dXdt{p{X,t)F{X,t)) = . 



We define an action of the gauge group on p with a plus sign instead of a minus sign when acting on 
F: 

Vip = {Di + Ai)p 
VtP = {dt + Q)p. 

This means that they p and F have opposite charges under the M gauge group, such that pF is neutral 
and Vt{pQ) = dtipQ). We have thus 



j dX {Vtp{X, t)F{X, t) + piX, t)VtF{X, t)) = 



Using equation Q for VtF{X, t) and integrating by part on the spatial directions, this equation 
becomes 

j dX (vtp{X, t) - ^V'VipiX, t) j F{X, t) = (6) 
This equation will be automatically satisfied if 

Vtp = ^V*V,p . (7) 

Moreover, if the market is complete equation Q must be true for all functions F{-,t) which imposes 
equation ([7|. This is the Kolmogorov forward equation, written in a covariant way. It should be 
noted that p{X, t) is a density, which means that the Levi-Civita connection does not reduce to a 
partial derivative as would be the case for a scalar. 



3 Asymptotic implied volatility 

We consider a stochastic volatility model where the variable is a forward price or rate F with a 
volatility variable V: 

dF = aF{F,V)dWi 

dV = fiv{V)dt + av{V)dW2 (8) 
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with (dVFidVF2) = pdt. 



Our computation of an asymptotic expansion at short time of imphed volatihty at strike K involves 
three steps: 

1. Compute an asymptotic value of the transition probability from initial state Fq, Vq at time 
to K, V at time t using a heat kernel expansion; 

2. Compute K^ap{K, V)5{F — K)~^ using a saddle point method; 

3. Integrate over time to compute the time value; 

4. Compare to the same formula for the Black-Scholes model to extract the implied volatility. 
3.1 Heat kernel expansion 

At short time the solution of equation ([T]) with initial condition p{X, 0) = 6{X — Xq) is asymptotically 
given by the heat kernel expansion |DeW65j (see |Vas03] for a review)]^ 

d'{Xo,X) 

piX,t) = -^^^A(Xo,X)P(Xo,X)e 2t ^ak{X,,X)t' (9) 

g{X) is the determinant of the metric at point X: 

9 = Det(5ij) . 

d{XQ,X) is the geodesic distance between the starting point Xq and the end point X, this is the 
minimal distance between Xq and X. It can also be written as 

d\Xo,X) ^^^^ f\tg,,X^X^ 

where the minimum is taken on all paths going from ^(0) = Xq to X{t) = X. (This is independent 
of t.) We denote by C this geodesic path. A{Xq,X) is the Van Vleck-Morette determinant 

1 W(Xo,X) \ 

2 dx^dxi J 




^g{Xo)g{X) 



V{Xq^X) is the parallel transport along the geodesic with respect to the M connection. It is such 
that its covariant derivative along the geodesic path is null: 



AiX'dt - / AidX' 
V{XQ,X)=e Jc =e Je 



Using Feynman path integral, the solution to equation ([7| can be written up to some normalization factor as 



p(X,t) oc j [DX]e 



At ( ]^g^jTP + A,X' 



where [DX] means integrating over all path X[s) going from X{Q) = Xo to X{t) — X. The normalization factor is 
the inverse of the same quantity with the integral computed over all paths with starting point Xo, so that the total 
probability is 1. It is generally not possible to compute this integral exactly. However it gives some hints on the 
asymptotic solution at short time: the solution will be dominated by the path corresponding to the minimal value of 
the integrand inside the exponential, which will be close to the geodesic path. 
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where the integral is computed on the geodesic path C. Finally, ai{XQ,X) are functions which are 
defined recursively with 

ao = 1 

and Oj's satisfy the differential equations 

Along a given geodesic curve parameterized by its geodesic distance from the origin d, this equation 
reads 

which can be integrated as 

fd 



d^ 

Note that when acting with the covariant derivative, the full charge under the M gauge group is 
carried by the parallel transport term V and that ^^^^'^^ is a scalar with respect to the Levi-Civita 

\/9{X) 

connection, a^'s are scalars with respect to both connections. 

In order to produce a first order expansion of the implied volatility, only the common multiplicative 
factor of expansion [9] is needed. In order to compute a second order term for the implied volatility, 
we wiU also make use of the first corrective term ait with 

3.2 Expected variance 

We now compute E [(T^(i^, V)6{F — ET)] . This quantity can be written as an integral over the terminal 
volatility variable V: 



¥.[al{K,V)5{F - K)] = j dV al{K,V)p{K,V; 



where p{F, V; t) is given by the heat kernel expansion Q with X = ^ ^ ^ and n = 2. The integrand 
can written as 

1 C-Dt + o(t) 

aUK,V)p{K,V;t) = —e t (11) 

with 

B = ^d{Fo,Vo;K,Vf (12) 

C = -2ln{aF{K,V))-^[ln{g{K,V)) + ln{A{Fo,Vo;K,V))]+A{K,V) (13) 
D = -ai{K,V) (14) 

where C is the geodesic curve joining (Fq, Vq) to (K, V), A is the integral of the M connection 

AiK,V) = -ln{ViK, V)) = j AidX' 
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and ai is given above as an integral over the geodesic path C in equation (10). 



The integral over (11) will be dominated at short time by the B term. More precisely, it will 
be dominated by the volatility Vmin which minimizes B{K, V) = ld{Fo, Vo; K, V)^. This is the final 
volatility which minimizes the distance between the initial conditions and the strike K. Expanding 
all functions in the neighborhood of Vmin, where -B'(Vtnin) = 0, the integrand can expanded as 



1 ---C-Dt --B"6V^ 
-e t e 2 



27rt 

1 /(^// _ ^,2\ 



lij(4) _ Ib^'^c) — + ^ + oit) + odd terms 

24 6 J t 72 t2 V ; 

where derivatives are with respect to V, where all functions B, C, D an their derivatives are taken at 
(K, V^in) and 5V = V^ — Vmin- When writing o(t), we have anticipated that after integration SV"^ ~ j. 
We have also anticipated that odd terms in 5V will not give contributions to the integral. 

Integrating over 5V , and using that the first even moments of the standard normal distribution 
are M2 = 1, M4 = 3 and Mq = 15, we get for the integral 



li?(4) _ 1^(3)^'' 

24 6 



3t 1 „|'o-,2 15t , , 
+ ^B^^^ — + o{t) 



B" 



72 



B"^ 



This can be rewritten as 



E[aUK, V)5{F-K)] 



R ^ ^ 

1 C-Dt + o{t) 



/27rt 



(15) 



with 



C 



D 



C+- HB") 



D + 



D + 



1 



2B" 
1 



C" - c"^ + 



15(4) ^(3) 5 ^^(3)- 



4 B" 



B" 



C" - c- 



15(4) 1/5(3)' 

+ ■ 



12 \ B" 
2" 



4 B" 



B" 



(16) 



(17) 



where all derivatives are with respect to V and all functions and their derivatives are taken at 

(if,^min). 



3.3 Time value 



The price of a Call of maturity T and strike K can be written as the payoff integrated against the 
risk-neutral distribution: 



Call(K, T) = e-''^ j dF{F- K)+p{F, T) . 



This can be written also as a double integral over strike and time as 



C8l\{K, T) = e 



{Fo-K)+ + JdF dt{F -K)+dtp{F,T) 



(18) 
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As is a forward, driftless, the Kolmogorov forward equation reduces to 

dtp{F,t) = ^d^F{aUF,t)p{F,t)) 
where ap{F,t) is the local (normal) volatility 

E[aUK,V)5{F - K)] 



(19) 



aUF,t)=E[aUF,V)\F = K] 



p{F,t) 



Plugging the Kolmogorov equation (19) in equation (18) and integrating twice by part on the F 

(20) 



variable, the Call price is finally obtained as an integral over time at strike K: 

rT 



Call(i^,T) = e-"^ {Fo-K)+ + \f Am[al{K,V)6{F - K)] 

L ^ JO 



1 



Using expression (15) for the integrand, the integral over time can be computed: 

T 







dm[al{K, V)5{F-K)\ 



t --r rr^ r ( [b\ D 

e t — y B erf c \ — 

vr I V t / 3 



^^{t-2B)e t + 2^3/2 erfcf'yi^ 



where erfc is the complementary error function, equal to the cumulative of the standard normal 
distribution up to \/2 factors: 

erfc(2;) = ^ J dye"?'' =2M(^-V2a 
The asymptotic expansion of this function at +oo 



erfc(3;) 



+ ^ + o 



X^/^T \ 2x2 Ax"^ 



with X = \/ Y gives the asymptotic expansion of the time value: 



„ / dtE[aUK,V)6{F - K)] = —= 

^ Jo ZV^TT 



B 



C-\n{B)-DT-—T + o{T) 



(21) 



3.4 Implied volatility 

The final step consists in computing the same expansion for the Black-Scholes model, which is simpler 
as there is no stochastic volatility to be integrated. The metric is given by the inverse of the variance: 

1 

9FF = ■ 

The Christoffel symbol is therefore 

^FF--J ■ 
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The M-connection components are computed using @ and ^: 

1 



Q 



2F 





1 


1 ^ 










a 


-TO 



The geodesic distance is 

d{Fo,K) 

The Van Vleck-Morette determinant is simply 

A(Fo,i^) = 1 . 

The parallel transport is 



V{Fo,K) = e 



K 



dFAi 



Putting all these elements together, the heat kernel expansion of p{K, t) is according to ([9]) 



p{K,t) 



1 



In 



2 K 



aKV2TTt\ K 
Multiplying by the local variance, we get 



2aH li-^t + o{t) 



Bbs 



E[al{K)6{F - K)] = a^K^p{K,t) 



Cbs - ^BS* + o{t) 



'2-Kt 



(22) 



with 



^BS 

Cbs 
5bs 



In^ 



K 



2a' Fo 
-ln{a)-^ln{KFo) 
a' 



(In fact formula ( |22[ ) is exact: there is no o{t) correction and it can be integrated exactly to get the 
Black-Scholes formula.) 



Writing equation (21) for both the stochastic volatility model and the Black-Scholes model, the 

(23) 



implied volatility is such that both quantities are equal: 



I + C + ln{B) + DT + ^T=^ + Cbs+ HBbs) + D^sT + T + o{T) 



Expanding the implied volatility a as a Taylor expansion 

a{K, T) = ao{K) + ai{K)T + a2{K)T^ + o{T^) 



and plugging this into equation (23) on the Black-Scholes side, we get 

2 \ 



2alT"' \FoJ 



In' 



+ ln 



(■ 


-2^T 


C2 

- 2 — 


+ 3 














1 




-2^ 


8 


T + 


2al 







] -ln{ao)-^T-lln{KFo) 

I (To 2 

^T=f + C + .n(B) + 5T+A 

-fo 
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Coefficients must be equal at each order in T, which gives our final expansion of the implied volatility. 
Power —1 gives the order implied volatility 



In 



/2B 



d{Fo,Vo;K,Vn,in) 

which was already obtained in |BBF02] and |HL05j . 

The first order correction is extracted from the constant term: 



(24) 



C + ln{aoVKFo) 



2B 



(25) 



Finally the 0(T) term gives the second order correction: 

2 



O"0 



3 fai 
2 Uo 



2B V 10 



(26) 



This gives our final result as the implied volatility expansion 



a 



i + ^r+^T2 + o(r2) 

O"0 ctq 



(27) 



We stress that this result is exact in strike: for a given strike, we have computed exactly the three first 
coefficients of the Taylor expansion. Moreover, contrary to other expansions, the order 1 expansion 
is extracted from the order expansion of the probability. This technique allows us to extract a 
second order term for the implied volatility from the order 1 term in the probability expansion. This 
method can be used to compute the Taylor expansion of implied volatility up to any order, although 
the computation becomes more complicated and involves integrals of increasing dimension: the 
coefficient of the heat kernel expansion involves k + 1 integrals. 



3.5 At the money 

The computation we have performed makes the implicit hypothesis that we are not exactly at the 
money: K ^ Fq. Otherwise, the dominant term in the exponential would vanish and we could not 
use the asymptotic expansion of the erfc function at infinity. Precisely at the money, we should use 
instead a Taylor expansion in 0. As the implied volatility surface is smooth, we just take the limit 
of formulas (24), (25) and (26) at — > i^o- If we perform instead the Taylor expansion of the erfc 



function at 0, we find only the two first orders 



o-o(i^o 

<7l 



-C{Fo) 



D{Fo) 



Careful Taylor expansions of all quantities at the money can be used to check that this is indeed the 



limit of equations (24) and (25). Moreover, it can be seen that these limit are conditions for formulas 



(25) and (26) to be convergent, as B goes to at the money (at order 2 in the geodesic distance. 



which means that the numerators must in fact vanish at order 2). 
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3.6 CEV volatility 



Instead of Black volatility, the asymptotic expansion can be computed for other local volatility models. 
Without stochastic volatility, the SABR model reduces to the CEV model. The local volatility part of 
the model is thus taken into account exactly without introducing approximation besides the stochastic 
corrections. In view of our application to the SABR model, we will compute here a CEV implied 
volatility. There are closed formulas for this model, involving Bessel functions. This implied volatility 
can therefore be used in the CEV pricing formula in order to get the price of the option. 

For a CEV model with parameter Pq and volatility factor a, such that 



the function B, C and D are 



dF = o-F*dVF 



Bo 
Co 

Do 



2a2 



In'(go) 



-ln{a) - -poHKFo) 
(3o{2 - (5o)a^ 



with 



° ■ /3o<l 



In 



l-/3o 
K 



/3o = l 



Formulas (24), (25) and (26) are modified as follows. 



l^ol 



/2S d{Fo,Vo]K,V^^) 



(28) 



02 



-70 



C + \n{ao) + UoHKFo) 



2B 



-) 




D + 2, 



/3o(2 - Po)al 



(29) 
(30) 



This gives the CEV implied volatility expansion 



0-1 



0"2 



1 + — T + — + o(T^ 
O"0 o-Q 



(31) 



The Black implied volatility formulas correspond to the special case f3o 
normal) implied volatility would correspond to /?o = 0. 



1. The Bachelier (i.e. 



3.7 Generalization 

This technique can be generalized easily to other parameterizations of the options prices. Consider 
a model with local volatility or stochastic volatility, for which there are closed form formulas for 
European option prices. It can be used as a proxy in the following way. 
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Denoting by Zi the parameters of the model, compute B^:{zi), C*(zj) and D^{zi), the quantities 
B, C and D of the asymptotic expansion (|21|) for this model at a given strike. 



Find parameters zf"^ such that = B (there can be several solutions). 

Choose a one-dimensional subset of the parameters Zj = Zi{\) which allows a wide range of 
option prices at the given strike and such that Zi{^) = zf^\ 

Compute derivatives of B^[zi)^ C^{zi) and D^{zi) with respect to A at A = 0. We use the 
notation B^. = S*(zi(0)), B'^ = dxB^X^iW) Ia=o • • • 

Write a Taylor expansion A(T) = AiT + X2T^ + o{T'^) and write the equality of the asymptotic 
expansion (21 ) for the model and the proxy model: 

B ~ , , ~ 3 B^ + XiTB'^ + XiT^Bi + lXiT^B'^ 
- + C + HB) + DT+—T=^ '-^ * ^ ' - 

+ C, + XiTCi + ln{B, + XiTB[) + D, + + o{T) . 

This gives the Taylor expansion of A: 

C-a+ln(B) - ln(5,) 



A 



BL 

D-D.-: 

A2 ■■ 



D — D^: — Xi^ — XiCl — jAf-B" 



B: 



• Plug parameters Zi{XiT + A2T^) into the closed form option price of the proxy model to get an 
approximate price of the option in the real model. 

The closer the models are, the better the approximation is. It is clear that if the proxy model is 
the real model itself, there are no corrections at all. This procedure consists in approximating only 
the differences between models at a given strike and not the option price itself. In the basic case of 
section |3.4| where the proxy model is the Black-Scholes model, the approximation leverages on the 
fact that the volatility surface is more regular than the option price. 



4 SABR Model 
4.1 Model 

The SABR Model |HKLW02] is a stochastic volatility model where the volatility is a local volatility 
function multiplied by a lognormal stochastic volatility: 

dF = VC{F)dWi 

dV = uVdW2 (32) 
with (dTyidW2) = pdt. The initial value for V is the parameteij^ a: 

a = V{0) . 

*We use the standard notation of a for the initial value of the volatility variable in the SABR model instead of Vq 
as in the previous section. 
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C{F) is a local volatility function, which is generally 



C{F) = 

(3 is a number between and 1 which controls the local skew. corresponds to a normal process and 
1 to a lognormal process. The implied volatility at time and at the money is the local volatility 
aFl^. 

For P<1, the origin F = can be reached with finite probability in finite time. If F corresponds 
to a positive variable, a boundary condition must be imposed. The asymptotic expansion does not 
distinguish between different boundary conditions, as the computation is local around the geodesic 
path. It is no longer valid when the probability of reflection or absorbtion at the origin modifies the 
probability distribution at the strike considered in a significant way. 

4.2 Order 0: metric 

In order to compute the order implied volatility, the only geometric object involved is the metric. 



According to the dictionary of section 2.2, its inverse is the covariance matrix 



V'C{Fy puV^C{F) 



This matrix is first simplified by changing the variable F to 

(33) 



^ dF 



C{F) 



which for C{F) = F^ reads for /3 / 1 



1-/3 

and for /? = 1 

In addition, we rescale the time such that u disappears of the equations while keeping the same 
solution of the equations (the variances which are the physical quantities are not changed) : 

t — > v'^t 
a 

a — > — 

V 

V — > 1 

At the end of the computation, the inverse transformation must be applied to the implied volatility: 

ov < — a 

The matrix in the set of variables {q, V) after this rescaling is 
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This is diagonalized by going from variables (g, V) to (x, y) with 

q-pV 

X = . 



The covariance matrix becomes 



and its inverse is the metric 



V . 



1/10 



9ij - .,2 I 1 



which corresponds to the infinitesimal distance 

ds 



2 dx^ + dy^ 

y2 



This geometry corresponds to the hyperbolic plane, in the Poincare half-plane representation [y > 0) 
|HLW05l IHLOSj . Geodesies are vertical lines and semi-circles orthogonal to the y = axis. The 
geodesic distance between two point {xi,yi) and {x2,y2) can be computed: 

dix,,yl; X,, y2) = cosh^ ( 1 + - ^^^^ + - V^)' 

In the {q, V) variables, going from q = 0,V = a to q,V the geodesic distance is 

.if q^ + (V -a)^ -2pq(V -a)\ 

For a given strike, i.e. a given q, it is minimized by the volatility 

Vrain = \/ ^ 2paq ^ q^ 

and the minimal distance is 



(i(0, a; g) = cosh ^{YlB^ — P " 



In 



{1 + p)a 



Equation (24) gives the order implied volatility 



-I) 



Klin + pa + q 



In 



{1+ p)a 

(We have dropped the absolute values as the numerator and the denominator have the same sign.) 

Plugging the expression for Vmin and going back to the original time, with v factors, the order 
implied volatility for the SABR model is 



^0 = . , ^ °^ ^ (34) 

/ V + 2pavq + v'^q^ + pa + qv 
m ' 



( 



K 



1 + p)a 
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with 



9 = S 



At the money, the hmit of this expression is simply 



which is the local volatility. 



4.3 Order 1: connection 

To compute the order 1 correction, we need the scalar factor in the time value expansion, given by 



C in equation (16), with C given in equation (13) 



For the hyperbolic plane, the Van Vleck-Morette determinant can be computed as a function of 
the geodesic distance: 

A - ^ 
sinh((i) 

We need also 

aF{K,V) = VC{K) = VK^ 



Using that d' = for the minimum distance at V = Vjoim we compute 

— — ^ 

aKiin(l - p2)sinh((i) 

Terms simplify against each other to give at y = Vmm 



C=-ln(y^^y~Kn+A (35) 



A is the integral of the connection 1-form on the geodesic. According to formula (|4]), the connection 
is given by 

A= ^ ^ , dln(C(F))- ^^^^-^I dV. (36) 
2(1 -p2) ^ ^ 2(1 -p2) ^ ' 



In fact, A can be rewritten in {x, y) variables as 



It must be integrated from 



yi J \ a 



to 



q-pV 



y2 J \ y 
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If /? = 1, the 1-form 

A = — , ^ dx 

is exact and can be integrated directly: 

We consider now the general case where /3 < 1. xi = X2, the geodesic is vertical line. As A is 
along dx, its integral is zero: 

^ = 



In other cases, the first part of equation (36) is an exact form and can be integrated directly: 

The second part must be integrated on the geodesic path. The geodesic is a semi-circle with origin 
(X, 0), radius R and going through (xi,X2) and (^1,2/2)- The origin is therefore 

xj-xj+y^- yl 
- 2{x,-x,) 

and the radius 



R = ^yl + {xi - X)2 . (39) 
We parameterize the geodesic by t = tan(^/2) where 6 is the angle on the circle: 

l-t2 

X = X + R 



l+t2 

y = ^iT^- 

In this parametrization, the geodesic distance is given by 

ds = dln(t) 

and we can compute 

- -A 41) - loit.) - (40) 



c2(l-p2) 2(1 \Fj (i_,3)^_p2 



with 

G{t) =tan~\t) 



V(a + 6X)2-(l-/3)2i?2 I ^(a + 6X)2-(l-/3)2i?2 I 



(a + bX)^ = (1 -/3)2i?2 



cii + t(a + 6(X -ii)) 



" + tanh-^ ;^R + tia + KX-R)) + 5^)2 < (1 - /3)2^2 



V(l - /3)2i?2 _ (a + 5X)2 yv^(l-/3)2i?2_(a + 5X)2^ 



(41) 
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a 



1-/3 







c = {l-(i)p 



and 



R-Xi + X 
R + Xi- X 



Summing equations (37) and (40), the integral of the connection is finally 

= [G{t2) - G{h)] 



Replacing A in equation ( 35 ) we get 



C = -^ln(aFo^y^inK^) + < 



(42) 



(43) 



p- 



(1-/3)^1^ 



[G{t2) - G{h)] (3 < 1 



(44) 



Restoring the factor u, the order 1 correction is given by equation (25): 



In' 



\J oP' ^ 2pauq + v'^q^ + pa + qv 
{1+ p)a 



(45) 



o"o is given by equation (34), C by equation (44) (with a divided by i^, also inside T^min)) where X, 
in 

/ 



i?, ti and t2 are given in equations (38), (39) and (42) from 

/ uq — p\J + 2pavq + v'^q^ \ 



yi 



\ 



—pa 
a 

V 



and 



X2 
V2 



\ 



\J c? ^ 2pavq + v'^q^ 



J 



and G{t) is defined in formula (41). 

Using this expression, the first order implied volatility is 

a = o-o ( 1 + — t + o{t) 

which is valid for all positive strikes. Exactly at the money, the formula we give must be replaced 
by its limit, which can be computed by a Taylor expansion or numerically. At the money and only 
at the money it appears to be equal to the original HKLW formula. This is not surprising as their 
expansion is in fact an expansion in both maturity and moneyness. 



01 



1 



2j^2(/3-l) , 1 _,,r,z;./3-l , 1 ,,2 1„2,,2 



+ -pavPF^'^ + -v' 



-p u 



(46) 
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4.4 Order 2 



To compute the second order correction to implied volatility, we need to compute D as defined in 
equation (17), with D = — oi defined in equation (10). 



We have to compute oi as defined in equation (10). Most of the integration can be done analyti- 



cally. We have first the integral of Q along the geodesic: 



AQ) 



Qds . 



According to equation (|5]), Q is 



Q 



1-/3 + 



2(1 - p2) ; F2(l-I3) 



Using the values defined in the previous section for X, R, ti, t2, a, b and c, its integral along the 
geodesic is 



13 + 



P 



with 
H{t) 



a + b{R + X)+ cRt 
(a + bX){l + t2) + bR{l - t^) + 2cRt 



R:^ H{t2) - H{h) 

[1 - /3)2/?2 _ (o + ln(t2) - Hh) 



(47) 



+ < 



cR 



V(a + 6X)2 - (1 - /3)2i?2 
cR 



V(l-/3)2i?2_(a + 6X)' 



tan 



: tanh 



cR + t{a + b{X - R)) \ 
y(a + 6X)2 - (1 - /3)2i?2 j 

.^f cR + t{a + b{X - R)) 
\^{l-l3YR'^-{a + bXf 



{a + bXf > (1 -/3)2/?2 
{a + bXf < (1 -/3)2i?2 . 



Note that in the denominator, the quantity ln(t2) — In(ii) is up to a sign the geodesic distance d. If 
/3 = 1, a^^^ reduces to 

(48) 



AQ) 



R \ X2 — Xi 



^ 8(1 -p2)^ • 

The Laplacian on the hyperbolic plane is in (x, y) coordinates 

D^D,=y\dl + d^y) . 



As the Van Vleck-Morette determinant A 



sinh(s) 



depends only on the geodesic distance s, its 



derivative on the orthogonal coordinate vanishes: d±A = 0. On the other hand, by definition the 
parallel transport on the geodesic curve has no covariant derivative along the curve: VgV = 0. As a 
consequence, there is no crossed term and both terms decouple: we have 

(the M charge is carried only hy V). 

The metric part can be integrated analytically |HLW05( IAvr07] : 



,(R) 



1 + 



1 / cosh((i) 1 



d \sinh((i) d 



(49) 



The last part to integrate is the M connection term V ^V^ViV. As the action of gauge transfor- 
mations on the heat kernel expansion is fully carried by the parallel transport term V, oi can only 
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depend on gauge-invariant quantities constructed from F = dA. We split therefore A = A^^^ + A^^^ 
into a pure gauge part and such that F = dA^^^: 



Aio) 



1 



dln(C(F)) 



2(1 



dln(C(F)) 



pC'jF) 
2(1 



dV . 



Forgetting A^*^^ which is pure gauge (it can be checked by hand that A^^^ will not contribute), we 
denote by V^^^ the A^^^ part of the parallel transport: 



p{i) = e"-^'^'=e-/c^« 



with 



(i-/3)7w 



[G{t2) - G{t,)] /5<1 



I 2(1 

Using A^^^ and A^^^ which in (x, y) variables is 



pin 



K 



(50) 



V + a] 13 = 1 . 



A, 



pC'jF) / p 

2 lyr^ 



dx — dy 



we can rewrite 



where the specific form of ^^^^ can be used to simplify terms 

a,4^) + dyA^'> = . 

Computing A^^^ and A^^') analytically, we compute numerically the x and y derivatives and inte- 
grate numerically along the geodesic curve to get 



2d 



dsy 



(51) 



This integral can be computed by a numerical quadrature with few points. For [3=1 the connection 
has in fact no curvature and therefore a\ =0. 

We compute also the following quantities (at V = Fmin): 

d" 



a{l- p2)l4^ij^sinh((i) 
B" = dd" 

m_ 3_ 

B" 

B" 



12 ^^„^cosh(d) 1 



(52) 



y2. 

mill 



sinh(d) d 



We need finally C and C" . Using our decomposition A = ^ In ^ -|- A'^^'^, we can write from formulas 
(pi) and dTel) 



C = --ln(i^F) 



-ln( 



d 



2 Vsh:ih((i) 



+ ^ln{B")+ln{l-p')+Ai . 



21 



Its first and second derivatives m V at V = l^nin are 

1^(3) 



2 B" 
1 



COS! 



h(d) 1\ 1^(4) 1 / 



2^ ivsinh((i) d) ^ 2 B" 2 \ B" 



We choose to differentiate A^^"^ numerically, by finite difference, as the analytical expression is very 
long and hard to simplify. 

Simplifying against — aj^^^ we get finally 



D = -a 



(Q) JA) 



1 1 



a) ' + - + 



aW" _ . — — 



1 ' 8 ' 2B" 

with af^^ given in equation (47), a"^'' in (51 ), B" in (52 ) and A^"^^ in (50) with G{t) defined in equation 



(41). For /? = 1, this expression can be simplified using equation (48) for a^'^'* and 



A^^y 







2(1 



. 



We have finally obtained the second order corrective term 



(72 _ 3 /CTI 
fJo 2 I (To 



^f5 + 3^-^ 
2B\ (To 8 



such that we get the quadratic approximation to implied volatility 

a = (Tofi + ^r+^r2 + (,(r2 

The computation has been done in redefined variables such that u = \. To restore the v factors, 
a must be replaced by ^, T by v^T and the final implied volatility must be multiplied by v. 

At the money, the formula for ^ looks divergent but its limit is well defined. We compute this 
limit numerically, although it could be done analytically. 



4.5 CEV volatility 



The results of section 3.6 can be used to invert he SABR volatility into a CEV fractional volatility. 



Using formulas of section 3.6 the implied CEV volatility is computed and used in the closed- form 



option prices of the CEV model. 

This proves to be useful for low (3 at low strikes or with very small volatility of volatility: only 
the corrections to the CEV model which come from the stochastic volatility are approximated, not 
the local volatility part. For example, at the money the first order coefficient of the Black volatility 
which is given by equation |46| becomes for the CEV volatility 



(i^o 



1 

— ; 
12 



1 



2 2 
P V 



The corrective term in o?T has disappeared. 
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4.6 Numerical results 

We present in figures [T] and |2] the implied volatility given by our expansion and compare it to the 
implied volatility computed by a two-dimensional finite difference method scheme. We also show 
for comparison the implied volatility given by the original formula of [HKLW02J. In this example, 
parameters are Fq = 4, a = 30%, f3 = .7, v = 40%, p = —.5 The FDM scheme is a second order 
Yanenko scheme with 400 points in strike, 200 points in volatility and 30 time steps. 

At very short maturities, all expansions are good. At first order our expansion is equal to the 
HKLW at the money but is more regular in strikes and is better in the wings as our computation 
does not involve any approximation in the moneyness. Our second order expansion is one order of 
magnitude more precis^ When maturity grows, first order expansions lose precision but the second 
order remain relatively good up to 10 years, where v^T = 1.6. At higher maturities, the second order 
expansion explodes quadratically and finally gives even negative volatilities at very long maturity and 
low strikes. At long maturities, a FDM or an other numerical method must be used, unless a valid 
long maturity expansion could be found. 

5 Conclusion 

We have presented a general method to compute a Taylor expansion in maturity of implied volatility 
of stochastic volatility models. We give exact formulas for the first and second order corrections. As 
an application, we have computed this expansion for the SABR model and compared it to the implied 
volatility given by a numerical scheme. It appears that it gives more precise results than the usual 
formula and extend the domain where a short maturity expansion can be used. Outside this range 
of validity, other methods must be used: numerical schemes or possibly other approximations. 

If a closer model with closed formulas than Black-Scholes exists, we provide a method to use this 
model as a proxy to extend the domain of validity of the expansion. It would be interesting to see 
the results of this method for the SABR model with a stochastic volatility model as a proxy. 

Obtaining exact option prices at all maturities would be a non-perturbative computation, which 
is a longstanding issue in theoretical physics. 
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^In fact at even shorter maturities, the FDM scheme we use is less precise than this second order expansion, especially 
in the wings where the probability density is very small. 
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2.5 years 



Orders - 
Order 1 - 
HKLW - 
FDM 





Figure 1: Implied volatility and relative error for the SABR model with parameters Fq = A, a = 30%, 
(3 = .7, v = 40%, p = —.5 and maturities 2.5 yr, 5 yr, 7.5 yr and 10 yr. On the left, implied volatility 
are plotted for our first order and second order expansions, the original formula of IHKLWO^ and 
are compared to the result of a FDM solution. On the right are the relative errors with respect to this 
reference solution. 
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15 years 



Orders - 
Order 1 - 
HKLW - 

FDM -- 



Order 3 - 
Order I - 
HKLW - 



20 years 




30 years 



Figure 2: Implied volatility and relative error for the SABR model with parameters Fq = A, a = 30%, 
(3 = .7, v = 40%, p = —.5 and maturities 10 yr, 20 yr and 30 yr. On the left, implied volatility are 
plotted for our first order and second order expansions, the original formula of IHKLWO^ and are 
compared to the result of a FDM solution. On the right are the relative errors with respect to this 
reference solution. 
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Appendix A Mean reversion 



At long maturity, the SABR model is not realistic as the volatility process is a geometric Brownian 
motion. In particular the variance of the volatility increases linearly in time. A direct extension 
would be to add mean reversion to the volatility process, either on the volatility or on the variance. 
The asymptotic expansion can still be computed at order 2. However its domain of validity is usually 
reduced: in addition to other conditions, the maturity must be small in front of the mean reversion 
characteristic time. 

We impose mean reversion on the volatility process as 

dV = uVdW2 + k{V - V)dt 

The metric is not modified as it describes the diffusion part. Expressions for A and Q are modified 

as follows: 

Pk{V-V) k{V-V) 

^ = ^SABR - ^^2(1 _ ^2)^/3 dF + ^2^2(1 _ ^2) 

1 K''{V-Vf 1 V 1 pI3k{V-V) 

Q = t^SABR + 7, 9\ + n'^ ~ '^l 



2 1/21/2(1 _ p2) ' 2 V 2v{l- p^)F^-f^ ' 
At first order, the integral of A along the geodesic is needed. Using the same notations as in 



section 4.3 where variables have been rescaled such that v = gets an additional term 

2(tan-Ht2) -tan-i(ti)) 



A A 

■A = AsABK + 



fV\ V V 
\a J V a 



The second order correction involves a one-dimensional integral which can be computed numeri- 
cally by a quadrature with a few points, although it may be possible to get an analytical expression. 
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